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Abstract 



A fruitful approach for solving signal deconvolution problems consists of resorting to a frame- 
based convex variational formulation. In this context, parallel proximal algorithms and related 
alternating direction methods of multipliers have become popular optimization techniques to 
approximate iteratively the desired solution. Until now, in most of these methods, either 
Lipschitz differentiability properties or tight frame representations were assumed. In this paper, 
it is shown that it is possible to relax these assumptions by considering a class of non necessarily 
tight frame representations, thus offering the possibility of addressing a broader class of signal 
restoration problems. In particular, it is possible to use non necessarily maximally decimated 
filter banks with perfect reconstruction, which are common tools in digital signal processing. 
The proposed approach allows us to solve both frame analysis and frame synthesis problems 
for various noise distributions. In our simulations, it is applied to the deconvolution of data 
corrupted with Poisson noise or Laplacian noise by using (non-tight) discrete dual-tree wavelet 
representations and filter bank structures. 
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1 Introduction 



Many works in signal/image processing are concerned with data restoration problems. For such 
problems, the original data y € £^(Z) is degraded by a stable convolutive operator T: £^(Z) — > 
£^(Z) and by a non- necessarily additive noise0 The resulting observation model can be written 
as z = T>a{Ty) where denotes the noise effect and a > is some related parameter (for 
example, a may represent the variance for Gaussian noise or the scaling parameter for Poisson 
noise). In this context, our objective is to recover a signal y G •^^(Z), the closest possible to y, 
from the observation vector z assumed to belong to £^(Z) and available prior information (sparsity, 
positivity,. . . ). In early works, this problem was solved, mainly for Gaussian noise, by using Wiener 
filtering, or equivalently quadratic regularization techniques. Later, multiresolution analyses were 
used for denoising by applying a thresholding to the generated coefficients [2]. Then, in order to 
improve the denoising performance, redundant frame representations were substituted for wavelet 
bases OH]. In [5l El 13 E] j authors considered convex optimization techniques to jointly address 
the effects of a noise and of a linear degradation within a convex variational framework. When the 
noise is Gaussian, the forward-backward (FB) algorithm [5] (also known as thresholded Landweber 
algorithm when the regularization term is an £i-norm [H [71 1^) and its extensions [HI [TO] can be 
employed in the context of wavelet basis decompositions and its use can be extended to arbitrary 
frame representations [11] . However, in the context of a non- additive noise such as a Poisson noise or 
a Laplace noise, FB algorithm is no longer applicable due to the non-Lipschitz differentiability of the 
data fidelity term. Other convex optimization algorithms must be employed such as the Douglas- 
Rachford (DR) algorithm |12) . the Parallel ProXimal Algorithm (PPXA) [13] or the Alternating 
Direction Method of Multipliers (ADMM) [141115]. These algorithms belong to the class of proximal 
algorithms and, for tractability issues, they often require to use tight frame representation for which 
closed forms of the involved proximity operators can be derived [13| [T^ [T5l [T6] . The goal of this 
paper is to propose a way to relax the tight frame requirement by considering an appropriate class 
of frame representations. 

In the following, we consider two general convex minimization problems, which are useful to 
solve frame-based restoration problems formulated under a Synthesis Form (SF) or an Analysis 
Form (AF). The SF can be expressed as: 



minimize fr{Lry) + ^Qs{x) (1) 

and the AF is: 



minimize (L^y) + 5s {Fy) ■ (2) 

F: £^(Z) £^(Z) (resp. F* : £^(Z) — )• £^(Z)) denotes the frame analysis (resp. synthesis) oper- 



^£^(Z) denotes the space of discrete-time real- valued signals defined on Z having a finite energy. 
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ator. For every r € {I, . . . , R}, fr'- — )• ] — cx3,+oo] is a convex, lower semicontinuous, and 
proper function, : e^{Z) is a stable convolutive operator, and for every s € {1, . . . , 5}, 

gs'. £^(Z) — 7> ]— oo,+oo] is a convex, lower semicontinuous, and proper function. In several works, 
SF has been preferred since AF appears to be more difficult to solve numerically \n\ \TE\ \T9\ [20] . 
In the proposed framework, both approaches have a similar complexity. 

This paper is organized as follows: in Section [2l the class of frames considered in this work 
is defined and their connections with filter bank structures is emphasized. In Section [3l we show 
how these (non necessarily tight) frames can be combined with parallel proximal algorithms in 
order to solve Problems ([T]) and ([2]). The proposed approach is also applicable to related aug- 
mented Lagrangian approaches. Finally, restoration results are provided in Section 3] for scenarios 
involving Poisson noise or Laplace noise by using Dual- Tree Transforms (DTT) and filter bank 
representations. 

Notation: Throughout this paper, To{'H) designates the class of lower semicontinuous convex 
functions if defined on a real Hilbert space 7i and taking their values in ]— oo,-|-oo], which are 
proper in the sense that their domain domc^ = {u G 7^ | (p{u) < +oo^ is nonempty. If (/9 G TqCH) 
has a unique minimizer, it is denoted by argmin(/7(n). The interior of a set C is denoted by int C. 

2 Frame representations 
2.1 Definitions 

Physical properties of the target signal y € i^C^), such as sparsity or spatial regularity, may be 
suitably expressed in terms of its coefficients x = {x{k))k£Z G i^C^) where y = Ylk^-oo'^i^)^k and 
{Gk)kez denotes a dictionary of signals in Such a dictionary constitutes a frame if there exist 

two constants /u and 7^ in ]0, +oo[ such that 

+ 00 

(yyee^iz)) iA\yf< \{y \ek)\^ <nyf- (3) 

fc=— oo 

The associated frame operator is the injective linear operator defined as 

{\/yG£\Z)) Fy={{y\ek)),^^, (4) 
the adjoint of which is the surjective linear operator given by 

+00 

{yx={x{k))^^^e£\Z)) F*x= x{k)ek. (5) 

fc=— 00 

When = F*, an orthonormal basis is obtained. Further constructions as well as a detailed 
account of frame theory in Hilbert spaces can be found in |21| . 
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A tight frame is such that, for some fi S ]0,+oo[, F*F = /il where I denotes the identity 
operator. In Condition ([3|), this means that fi = JI = /i. A simple example of a tight frame is the 
union of // orthonormal bases. Other examples of tight frames can be found in [22\ [23l 124) . 



2.2 A class of non necessarily tight frames 



We consider a linear operator F which is basically obtained by cascading a non necessarily max- 
imally decimated filter bank and a semi-orthogonal transform. A linear operator U: 

(£^(Z))'^, with N £ W and Q G N*, is said to be semi-orthogonal if there exists G ]0, -|-oo[ 
such that U*U = nij I . Recall that an analysis filter bank can be put under its polyphase form by 
performing a polyphase decomposition II/) followed by a real MIMO (Multi-Input Multi-Output) 
fihering V [25|: 

• the polyphase decomposition 11/) is an operator from to (^^(Z)) with D such that, 
for every y = {y{n))^^^ G -^^(Z), ^dV = {y^^^)i<j<D where y(^) = {y{Dn - j + l))^^^ is the 
j-t\i polyphase component of order D of the signal y. The adjoint operator of II/) is given by 
U}) : (£2(Z))^ £'^{Z) : {y^^'^)i<j<D ^u = {u{n))^^^ where, for every n G Z and j G {1, . . . , D}, 
u{Dn — j + 1) = y^^\n). So, 11^ allows us to concatenate D square summable sequences into a 
single one. It can be noticed that IIJjII^) = I and XIj^II^ = I , which means that YIq is an isometry 
and U-^^ = n^. 

• The MIMO filter V is defined as 



V 



'Vi. 



Vi. 



D 



(6) 



Vn. 



D 



where, for every i G {1, . . . , iV} and j G {1, . . . , D}, Vij : is a SISO (Single-Input 

Single-Output) stable filter. Hence, the impulse response of this filter belongs to ^^(Z) and its 
frequency response Vi,j is a continuous function. In addition, it is assumed that V is left invertible, 
that is: for every v G [—1/2, 1/2], the rank of the matrix v{i/) = \vij{i')]Ki<:N,i<j<D is equal to 
D. The adjoint operator of V is the D x N MIMO filter given by 



V* 



V* 

^1,1 



V* 



V, 



N,l 



V, 



N.D. 



(7) 



where, for every i G {1, . . . , A^} and j G {1, . . . , D}, V*j is the SISO filter with complex conjugate 
frequency response v* j. 

We have then the following result (the proof is provided in Appendix 15. ip : 

Proposition 2.1 The operator F = UqUVHd is a frame operator with frame constants fi = 
inf^g[_i/2,i/2] o"i„in('^) andfl = sup^g[_i/2,i/2] o"max(i^), where, for every v, arain{v) G ]0,+oo[ and 
'7max('^) S ]0, -|-cxd[ are the minimum and maximum eigenvalues of In addition, we 

have: 

F*F = fxun*DV*vnD. (8) 



^The notation v{iy)^ corresponds to the transconjugate of v{i 



4 



The resulting frame is not necessarily tight. However, when (Vz^ € [—1/2, 1/2]) v{i') v{u) = lu , 
the frame is tight. This includes paraunitary systems V as particular cases when D = N . 

Below, we provide examples of popular transforms in signal processing which belong to the class 
of considered frame representations. 



2.2.1 Example 1 — Dual-tree transforms 



In order to obtain low redundancy representations, frames such as the 2D M-band DTT have 
been proposed |26| [271 128) . The real (resp. complex) DTT consists of performing N = 2 (resp. 

= 4) M-band orthonormal wavelet decompositions in parallel where, for every i € {1, . . . , A^}, 
J7j : £^(Z) — )■ £^(Z) denotes the i-th orthogonal wavelet transform. An orthogonal combination of 
the subbands modeled by (£^(Z)) (^^(Z)) is applied to ensure directionality properties. 
{Ui)i<i<N and $ are related to U by the following relation: 



C/ = $ 



Ui 

C/2 









Un 



(9) 



Since <I)*<I) = = I and, for every i G {l,...,iV}, U*Ui = UiU* =1, ?7 is an orthogonal 

transform. In the general form of DTT, each orthonormal wavelet decomposition is preceded by a 
prefiltering stage related to the discretization process. More precisely, it takes the form: 



(10) 



The variables Q and D as defined before are thus equal to N and 1, respectively. The left 
invertibility condition here reduces to the fact that the frequency response + • • • + |i^Ar,ip 

does not vanish. Note that, due to the presence of prefilters, the discrete DTT is not a tight 
frame in general. Furthermore, an extended class of DTTs can be obtained by using tight frame 
overcomplete representations (f/j)i<j<7v having the same frame constant ^jj (e.g. redundant wavelet 
representations derived from orthogonal filters after appropriate renormalization), in which case U 
as defined by ([9]) is a semi-orthogonal transform. 



2.2.2 Example 2 Filter banks 

Analysis filter banks with perfect reconstruction correspond to the case when U = 1 and Q = N . 
D is the decimation factor and N is the number of channels. The redundancy introduced by such 
a filter bank structure is N/D > 1. More details about filter bank design can be found in |29| . 
Lapped transforms |30) can also be implemented with filter bank structures. As already mentioned, 
filter banks do not constitute tight frames in general. Note that if D = 1, a fully undecimated filter 
bank is obtained. 

Tight frame representations have been widely used for data recovery by using convex optimiza- 
tion methods. In the next section, we recall some convex optimization tools and show the relevance 
of the considered class of frames in recent optimization approaches. 
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3 Use of parallel proximal algorithms 



A number of algorithms such as the alternating split Bregman algorithm |31| . augmented La- 
grangian techniques |14] and parallel proximal methods |32) have been recently proposed to ad- 
dress possibly nonsmooth convex optimization problems encountered in the solution of restoration 
problems. Here, we will focus on this specific class of algorithms which have proven to be useful 
for solving problems like ([1]) and ([2]) when a tight frame representation is employed. A common 
feature of the aforementioned algorithms is that they require a large-size linear system to be solved 
at each iteration. When non tight frame representations are used, the computational cost of the 
associated inversion may become prohibitive. We will see however that for the class of frames 
introduced in Section [21 this inversion can be performed in an efficient manner. This fact will be 
demonstrated by considering an extension of the Parallel ProXimal Algorithm (PPXA), hereafter 
designated as PPXA+. As shown in ^32j . PPXA+ is closely related to the alternating direction 
method of multipliers and its parallel extensions |15) . 

We recall that the proximity operator [33j of a function ip ^ TQ{7i) is defined as 

1 2 

prox : 7^ ^ : ?; 1-^ argmin — ||u — w|[ + ip{u). (11) 

ueH 2 

It can be observed that when if is the indicator function lq of a nonempty closed convex subset 
C of H, i.e. it takes on the value in C and -|-oo in 1-L\C, prox^^ reduces to the projector Pc 
onto C. Other examples of proximity operators corresponding to potential functions of standard 
log-concave univariate probability densities have been listed in [5l|TTl|T3]. The proximity operators 
employed in the experimental part of this paper, are recalled below. 

Example 3.1 (soft-thresholding rule) 

Let a > 0, and set cp: M — )• M: ^ i-^ Oi\^\. Then, for every E M, prox^^ = sign(^) max{|^| — a, 0}. 

Example 3.2 (Poisson minus log-likelihood function) |11) 
Let a > 0, X ^ 0, and set 



Then, for every ^ € 




"C; if X> and^ > 0; 

if X^O and^>0; (12) 
otherwise. 



prox .= ^-°+^^° £±i^. (13) 



One of the difficulties in the resolution of problems like ([TJ or ([2]) is that they involve linear 
operators. Unfortunately, the proximity operator of the composition of a linear operator and a 
convex function takes a closed form expression only under restrictive assumptions as stated below: 
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Proposition 3.3 [E] 

Let Q be a real Hilhert space, ip € ^o{Q) (^nd let L: % Q denote a hounded linear opera- 
tor. Suppose that LL* = for some x £ ]0,+oo[. Then, 99 o L € ro(^) and prox^^^ = 
I +x-'L*ipvox^^-l)L. 

When dealing with hnear operators L such that LL* 7^ ^or any x > 0, proximal algorithms 
requiring the inversion of a linear operator at each iteration can however be designed. For ex- 
ample, Algorithm [1] (resp. Algorithm [2]) can be applied to Problem ([1]) (resp. Problem ([2])). (In 
these algorithms, the sequences {ar/)i<r<R and (6s/)i<s<s' model possible numerical errors in the 
computation of the proximity operators at iteration i.) 

Algorithm 1 

iiitialization 

{Vrh<r<R G ]0, +00[^ , (k«)i<s<5 G ]0, +Oo[^ ; Ko)l<r<i? G (^^(Z))^, K,o)l<s<5 G 

XQ = argmin„g^2(^) Ylr=i rir\\LrF*u - Vrfi\\ + Y,s=i ''^s\\u - Wsfi\\ 

For ^ = 0,1,... 

For r = 1, . . . , pr/ = proxjj^^ Vr/ + Or/ 

For s = 1, . . . , 5 rs/ = prox^^/^^ w.,/ + hs/ 

A, G]0,2[ 

V^-R II T t:^* 1(2 , II ||2 

ci> = axgTaia^^f2^^-^}_^^^^r]r\\LrF u - pr/\\ + 2^s=i'^s\\u - rs/\\ 
For r = 1, . . . , Vr,i+i = Vr,e + Xe{LrF*{2ce - xe) - pr/) 

For s = 1,. . . ,S Ws/+i = Ws,e + Xi{2c(, - xe - rg/) 

Xi+i = Xi + \i{ci - Xg) 



Algorithm 2 

. nitialization 

(??r)l<r<R G ]0, +Oo[^ , (Ks)i<s<5 G ]0, +Oo['^ ; {Vr,o)l<r<R G (^^(Z)) , (Ws,o)l<s<5 G (^^(Z)) 

yo = argmin^g^2(g) 'l2r=i^r\\LrU - Vrflf + Ef=i - Ws,of 

For £ = 0, 1, . . . 

For r = 1,...,R p^.f = proxy^/^^ w^^f + 

For s = 1, . . . , 5 rs,f = prox^^/^^ w^^f + 65,^ 

G]0,2[ 

Q = argmin„g^2(2) Ylr=l Vr\\LrU - Pr/f + J2s=l'^s\\Fu - Ts/f 

For r = l,...,R Vr/+i = Vr,e + \i{Lr{2ci - ye) - Pr/) 

For s = 1, . . . , 5 Ws,£+i = + A<?(F(2q - y^) - rg/) 

yi+i = ye. + Xe.{ci> - Vi) 



The convergence of the sequence (x^)^gz (resp. (y^)^gg) generated by Algorithm [T] (resp. Al- 
gorithm [2]) to an optimal solution of the related optimization problem is guaranteed under the 
following technical assumptions (see |32| for more details): 

Assumption 3.4 

i. ( n^Li int dom fr o LrF*) n ( flfj/ int dom Qs) fl dom 55/0 
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(resp. ^ int dom o j PI ^ P|f^]^ int domgis o Fj Pi dom (75 o F 7^ 0j. 

2. There exists A g]0, 2[ such that € N) A < A^+i < Xi>. 

3. (Vr G {1, . . . , R}) llflr.dl < +00 and (Vs G {1, . . . , 5}) Et°S < 

Note that quadratic minimizations need to be performed in the initiahzation step and in the 
computation of the intermediate variable q at iteration i. These amount to inverting a large-size 
linear operator which is untractable for arbitrary frames. We will now show that the class of frame 
operators F introduced in Section [2.21 allows us to overcome this difficulty. 

To do so, recall that the convolutive operators (Lr)i<r<i? can be put under a polyphase form 

m 

(Vr G Lr = W,IlD (14) 

where the polyphase decomposition operator Ylo has been defined in Section 12.21 and Wr = 
[Wr,i, Wr,D] is a MISO (Multi-Input Single-Output) filter (for every j G {1, . . . , D}, Wrj : 
£^(Z) is a SISO filter). Now, invoking Proposition 12.11 and making use of Sherman-Morrison- 
Woodbury identity yield: 
• for SF (Algorithm [1]) 

R R R 

^7]rFL;LrF* + Kiy = I -F ( ^ r/^L^L^) I +F*F ( ^ T^^L^L^) ) ~ F* 

r=l r=l r=l 

R ^1 
= 1 -FIil{Y.'^rW;Wr)TiD{t^l +^JiumDy*y[Y.'^rW;Wr)JlDy F* 
r=l r=l 
R R 
= 1 -FUh{^VrW;Wr)Ul +f,uV*V{^r]rW:Wr)) H^F*, 



r=l r=l 



(15) 



• for AF (Algorithm [2]) 

R R 

(^^r]rL;Lr + KF*Fy =U*o(^Y,^rW;Wr + KfiuV*Vy Un, (16) 

r=l r=l 

where k = X^f=i^s- Note that the idea of using the Woodbury matrix identity to handle the 
inversion for SF was proposed in [34] in the context of tight frames. The inversions in (jlSp (resp. 
(jl6p ) can be performed by noticing that, as {Wr)i<r<R and V are multivariate filters with frequency 
responses {Wr)i<r<R and v, k1 +fJ^uV*V {^^^^rirW*Wr) (resp. Ylr=i^rW*Wr + k^c/^*^) is a 
MIMO filter with frequency response: for every v G [—1/2, 1/2], 

R R 

r=X r=l 

Hence, by resorting to Fast Fourier Transform implementations, the problem reduces to the inver- 
sion of D X D matrices. 



8 



4 Experimental results 



We apply the proposed optimization method to the restoration of images degraded by a blur and 
a Poisson (resp. Laplace) noise. For simplicity sake, the proposed formulation has been presented 
in the ID case, but its 2D extension is straightforward. The scaling factor of the Poisson (resp. 
Laplace) noise is a = 0.1 (resp. a = 10). The data fidelity term is related to the Poisson (resp. 
Laplace) minus log-likelihood, which corresponds to the potential function ip defined by (|12p (resp. 
by Example [XT])- Li the considered examples, the regularization term simply is an £i-norm modeling 
the sparsity of y through a frame representation. Two different frames are considered: a (Complex) 
Dual- Tree Transform, presented in Section [2.2.11 and an eigenfilter bank decomposition |35| I36|. 
which is a special case of the subband structure presented in Section 12.2.21 We constrain the data 
values to belong to [0,255], so defining a closed convex constraint set C. The considered SF (resp. 
AF) problem is a particular case of Problem ([T]) (resp. ([2])) where R = 2, S = 1, fi = (p, f2 = ic^ 
and gi = t\\ ■ ||i with r > 0. The last function corresponds to the regularization term operating 
in the frame domain. In the considered problems, Li = T and L2 = I • The proximity operators 
associated to /i, /2, and gi are derived from Example 13.21 the projection onto C, and Example 13. 11 
In our simulations, the parameter r is empirically chosen to maximize the signal-to-noise-ratio 
(SNR). In general, the value of this parameter is not the same for FA and FS. 

Note that an alternative to the proposed approach consists of resorting to primal-dual algorithms 
|37l [38] l39l l40j |4T| . These algorithms are appealing as they do not require any operator inversion 
and they can thus be employed with arbitrary frames. However, after appropriate choices for the 
weights {rir)i<r<R, (ks)i<s<5, and the relaxation parameter Xi, the proposed approach appeared 
to be faster. For example, similar frequency domain implementations of the Monotone + Skew 
Forward Backward Forward (M+SFBF) algorithm |41) were observed to be about twenty times 
slower than the proposed method in term of iterations and computation time. 

Figure [T] shows the restoration results for a cropped version of the "Barbara" image in the 
presence of Poisson noise and a uniform blur of size 5x5. We adopt a SF criterion and we consider 
a tight version (i.e., for every i £ {1, ... ,4}, T^.i = I) as well as a non-tight version of Complex 
DTT. The Complex DTT p8| is computed using symlets of length 6 over 3 resolution levels. In 
order to efficiently perform the inversions in (fTSl) and (I16p . fast discrete Fourier diagonalization 
techniques have been employed. The use of the non-tight Complex DTT including prefilters allows 
us to improve the quality of the results both visually and in terms of SNR and SSIM |42| . Figure [2] 
displays a second restoration example for a cropped version of the "Straw" image in the presence 
of Laplace noise and a uniform blur of size 5x5. AF results are presented by using a DTT and an 
eigenfilter bank {D = 8 and N = 14) computed from the degraded image. This formulation leads 
to better results than those obtained with SF. Significant gains in favour of the eigenfilter bank 
can be observed . 
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Original Degraded Tight-frame Complex DTT Complex DTT 

SNR = 11.4 dB SNR = 13.3 dB SNR = 14.2 dB 

SSIM = 0.53 SSIM = 0.69 SSIM = 0.73 



Figure 1: Cropped versions of Barbara image (size 128 x 128). Restored images using SF and 
complex DTT. 




Original Degraded DTT Eigenfilter bank 

SNR = 14.8 dB SNR = 16.7 dB SNR = 17.3 dB 

SSIM 0.42 SSIM = 0.64 SSIM = 0.69 



Figure 2: Cropped versions of Straw image (size 128 x 128). Restored images using AF considering 
DTT and eigenfilter banks. 



10 



5 Appendix 



5.1 Proof of Proposition 12.11 

Eq. ([5]) is a direct consequence of the fact that U is semi-orthogonal and I\*q is an isometry. By using 
Parseval's formula, for every y € -^^(Z), 

llFyf = tiu\\V\lDV\? = f'^ \\v{v)y{y)\\^dy (18) 

where y{v) ~ [y'-^'(i^), . . . ,^'■■^■'(1')]^, and, for every j <E {1,...,!?}, y^^^ is the Fourier transform of the 

discrete signal y^^K By using the fact that, for every G [—1/2,1/2], o'min(j^)||y(i^)|P < ^ 
we get 



.1/2 .1/2 

^iu inf ^a^W) \\yily)fdly<\\Fy\\^<^lu sup a^^y) Wvi^Wdv. (19) 

i^e[-i/2,i/2] J-1/2 i>e[-i/2,i/2] J-1/2 

Besides, we have 

»l/2 -P ^1/2 +00 +CX3 



"1/2 7 = 1-^-1/2 7=0 n=-oo n=-oo 



2/11 ■ (20) 



Now, it can be noticed that since v is continuous, fJmin and fJmax are continuous functions too. 
By invoking Weierstrass theorem, there thus exist 1/ G [—1/2,1/2] and F G [—1/2,1/2] such that 
inf^g[_i/24/2] crmin(i^) = cr^in{v) and sup^,g[_i/2,i/2] crmax(i') = o-max(t^)- Finally, (Tinin(i^) > since 
rank {v{i/)^v{ii)^ — rank (■w(ii)) = 
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